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In the most prominent current scenario of galaxy formation, galaxies form hierarchi- 
cally through the merger of smaller systems. Such mergers could leave behind dynamical 
signatures which may linger long after the event. In particular, the globular cluster sys- 
tem (GCS) of a merging satellite galaxy may remain as a distinct sub-population within 
the GCS of a massive galaxy. Using the latest available globular cluster velocities and 
■ metallicities, we present the results of a search for grouping in the GCS of our nearest 

large spiral galaxy neighbor, M31. A modified friends-of-friends algorithm is used to 
identify a number of possible merger remnants in projected position, radial velocity and 
[Fe/H] parameter space. Numerical simulations are used to check that such merger rem- 
nants are indeed plausible over the timescales of interest. The identification of stellar 
streams associated with these groups is required in order to confirm that they represent 
merger remnants 

> 

1. Introduction 

In hierarchical structure formation scenarios, galaxies form at least in part from the merger 
and accretion of smaller clumps of matter (e.g., Searle & Zinn 1978; Cole et al. 2000). Clumping in 
the early universe is believed to have resulted from small perturbations in the primordial density 
field, instabilities that triggered the growth of structure from the "bottom up" . Such primordial 
clumps would then evolve into dwarf galaxies like those observed in the present epoch, or may 
merge to form more massive galaxies. 

There are indications that such mergers are ongoing, as is evident from observations of the 
Sagittarius dwarf galaxy's current interaction with the Milky Way (Ibata, Gilmore & Irwin 1995; 
Dohm-Palmer et al. 2001; Newberg et al. 2002). There are also indications of tidal distortions 
in M31's dwarf satellite NGC 205 (Walterbos & Kennicutt 1987; Bender, Paquet & Nieto 1991). 
Furthermore, a wide-field imaging survey of M31 by Ibata et al. (2001) has revealed the presence of 
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a giant stream of metal-rich stars within the halo of M31. This feature may have originated from 
tidal interactions with M31's dwarf companions, NGC 205 and M32. Either the process of satellite 
accretion is not uncommon, or we are observing our Galaxy and M31 at particularly unique stages 
of their evolution. Recent studies by Reitzel & Guhathakurta (2002) and Ferguson et al. (2002) 
have also found signs of substructure and tidal trails in the M31 halo which support the possibility 
of past accretion events. 

Globular clusters (GCs) have been identified within some of the dwarf satellites of the Milky 
Way (Da Costa & Armandroff 1995; Hodge et al. 1999; Oh, Lin & Richer 2000) and M31 (Grebel, 
Dolphin & Guhathakurta 2000), as well as in the dwarf members of distant galaxy clusters (e.g. 
Durrell et al. 1996; Miller et al. 1998). Based on their positions, distances and radial velocities, there 
are four Milky Way globular clusters which are believed to be associated with the Sagittarius dwarf 
(Ibata, Gilmore & Irwin 1994; Ibata et al. 1997; Bellazzini, Ferraro, & Ibata 2003). Ultimately, 
these clusters — M54, Arp 2, Terzan 7 and Terzan 8 — will have their orbits randomized by the 
strong potential of the Milky Way and they will presumably take their place as bona fide members 
of the Galactic globular cluster system. 

Therefore, in the context of hierarchical galaxy formation, one would expect that the globular 
cluster systems (GCSs) of massive galaxies could include the cannibalized globular clusters of past- 
accreted dwarfs. For instance, this seems to be a plausible mechanism for building up the GCSs of 
giant elliptical galaxies (e.g. Cote, Marzke, & West 1998). Such a scenario raises the question of 
whether the globular cluster systems of merged dwarfs would remain distinct dynamical entities or 
instead become hopelessly mixed with the larger galaxy's GCS. If they remain recognizable over 
significant time-scales, satellite GCS remnants could provide us with a powerful dynamical probe 
of the formation and evolution of galaxies. 

The first significant effort at distinguishing sub-groups within the GCS of an external galaxy 
was performed for the M31 system by Ashman & Bird (1993). These authors suggested that 
an apparent sub-clustering in the positions and velocities of M31 globular clusters represents a 
surviving signature of progenitor gas clouds from which the galactic halo itself may have formed. 
Ashman & Bird employed a technique in which globular clusters were grouped based on deviations 
from global mean velocities and the velocity dispersions between each cluster and its N nearest 
neighbors. Their sample included 144 M31 GCs with spectroscopic data from Huchra, Brodie & 
Kent (1991), many of which had large velocity uncertainties (> ±50 kms -1 ). 

In many respects, the M31 globular cluster system provides an ideal target in the search for 
substructure and grouping. More than 435 confirmed globular clusters have been identified in 
M31 (Barmby et al. 2000). A recent spectroscopic study of the M31 GCS by Perrett et al. (2002) 
has now pushed the number of clusters with radial velocities and metallicities beyond 300. In 
addition, Perrett et al. (2002) quote typical velocity uncertainties for the majority of the clusters 
of ±12 kms -1 , a significant improvement over many previous studies. The time is right to take a 
closer look for signatures of hierarchical formation within the M31 globular cluster system. 



-3- 



There are several important things to consider in this type of study. One must define what is to 
be regarded a significant "grouping" in the available parameter space. In other words, appropriate 
grouping criteria must be established a priori. Furthermore, one must investigate the likelihood 
that any observed groups may be mere coincidental associations rather than actual substructure. 
A greater difficulty is the interpretation of the significance of the grouping results: do the groups 
represent true merger remnants that have survived to the present epoch? 

In this paper, we develop a technique to identify potential dynamical remnants within a galaxy's 
globular cluster system. This group-finding algorithm is described in Section 2. In Section 3, we 
apply the group-finding technique to the sample of M31 GCS positions, velocities and metallicities 
from Perrett et al. (2002). We identify several potential groups of clusters in this sample. Monte 
Carlo simulations are performed on randomized globular cluster populations to characterize the 
significance of the identified groups. Section 4 describes our use of N-body simulations to investigate 
whether these groups could represent actual merger remnants. Finally, the results of this study are 
summarized in Section 5. 

2. The Technique of Group Identification 

In most applications, simulated data has the advantage over observational data in that it 
provides full three-dimensional positional and dynamical information. Observational data is, of 
course, less ideal: we must cope with projected positions and line-of-sight velocities and we can, at 
best, make generalized assumptions about the true dynamical behavior of any system. Additionally, 
every observed parameter carries with it an uncertainty, the value for which may vary from object 
to object even within the same dataset. These are some of the considerations that must be made 
when devising a group-finding algorithm for observational data. 

2.1. Group-Finding Algorithm 

The basic technique used in our group identification is a friends-of-friends algorithm. This 
algorithm links together particles that fall within a specified distance in parameter space. Each 
distinct set of joined particles constitutes a group. In the case of the M31 globular cluster system 
observations, the parameter space is defined by X and Y (projected) positions, radial velocities v r , 
and metallicity [Fe/H]. 

The appropriate weighting of these parameters must be determined in advance in order to set 
the group-linking criteria. Since each parameter has units, each must be converted to an equivalent 
unitless quantity on the same scale before calculating distances in parameter space. For example, if 
X and Y each span a range of 100 arcmin, the velocities range over 600 km s _1 , and metallicity spans 
only 3 dex, an appropriate metric must be devised such that the group is not defined principally 
based on [Fe/H] owing to its comparatively small numerical range. In essence, the metric puts all of 
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the parameters on to the same effective scale. In addition to the metric, the group-finding algorithm 
requires a linking length which defines the maximum allowable distance between successive group 
members in the specified parameter space. 

Therefore, for each GC, the algorithm links a nearby GC to it if 

(A S ) 2 = J>|(A^) 2 <£ 2 , (1) 

i 

where As is the total separation (in parameter space) of the neighbor to the GC in question, ji is 
the scaling or weighting factor defined for parameter i, Axj is the separation between the GC and 
its neighbor in parameter i, and L is the linking length. Traditionally, C is specified as a fraction, 
A, of the mean inter-particle spacing for the entire system: 

C = AAs (2) 

The advantage of using a friends-of-friends technique over average density methods is that it can 
follow elongated structures, features that may be expected in tidal disruptions. 

The error bars associated with each input parameter are accommodated by the addition of an 
"uncertainty ellipsoid" surrounding each GC in the available parameter space. A GC's parameter- 
space neighbor is linked with it if the minimum distance between the uncertainty ellipsoids is less 
than the specified linking length. 



2.2. Metric and Linking-Length Determinations 

In order to define a suitable metric to be used in the definition of a group, a method was 
devised to test the group-finding algorithm on simulated associations using a range of scalings 
and linking lengths. Groups of objects with various population sizes, spatial distributions (e.g. 
spherical, linear and curved streams), galactocentric positions and orientations were created with 
general parameters similar to those of Local Group dwarf galaxies. The simulated groups were 
given internal velocity dispersions of a v = 6 to 12 kms -1 , mean metallicities of ([Fe/H]) = —1.2 to 
-1.6 dex, and [Fe/H] spreads of cr [Fc/H] = 0.3 to 0.5 dex (Mateo 1998; Da Costa et al. 2000). 

The GCs in these simulated groups were then superimposed over a background designed to 
mimic a realistic distribution of M31 globular clusters. This background of several hundred GCs 
consisted of the following components: 

1. a random, metal-poor, spherical halo distribution, 

2. a rotating, metal-rich thick disk distribution inclined at 12.3° to the line-of-sight, and 

3. a rotating, intermediate-metallicity ellipsoidal bulge distribution. 
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Projected positions, radial velocities and metallicities for the fabricated group members and 
additional background objects were analysed using the group-finding algorithm. The procedure 
was cycled through a range of linking lengths and parameter scalings to determine the optimal 
combination for group retrieval. A minimum group membership of uqq = 4 clusters was specified 
to define a group. 

Setting the X and Y scaling to 7x = 7y = 1-0 arcmin^ 1 to define the baseline, it was found 
that the best recovery rate occurred when A = 0.1 and the velocities were scaled by 7^ = 0.13 
km -1 s for cases when metallicities were not considered in the group recovery. With the addition of 
metallicity information, the optimal velocity and [Fe/H] scalings were found to be j v = 0.14 km -1 s 
and 7[Fe/H] = 9-7 dex -1 , respectively, with a linking length given by A = 0.09. These optimal linking 
length and parameter scaling values were used to search for grouping in the globular cluster system 
of M31. 

3. Substructure in the M31 GCS 

3.1. The GCS Sample 

The best available spectroscopic sample of > 300 GC positions, radial velocities and metal- 
licities as defined in Perrett et al. (2002) was examined for signs of grouping. Quoted errors on 
individual velocities and metallicities were used to define the uncertainty ellipsoids about each data 
point in parameter space. Uncertainties on the target positions were typically around 0.2 - 1.0 arc- 
seconds; these were deemed to be negligible and thus position uncertainties were not incorporated 
in the analysis. 

The group-finding algorithm was used to analyze the data based on the parameter scalings and 
linking lengths determined in Section 2.2. Again, a minimum of four clusters was required in order 
to define a group. A search for groups was conducted in position/velocity/metallicity parameter 
space for the sample of 301 globular clusters with published spectroscopic metallicities. Roughly 
one dozen groups were initially flagged within the GC population. The robustness of these groups 
was then tested against small (~ 10%) changes in the specified parameter scalings and linking 
length. Groups that were not found to be stable against these perturbations were rejected. 

At this point it should be noted that the inclusion of metallicity information may not necessarily 
improve the effectiveness of the group-finding results. The [Fe/H] uncertainties for many of the M31 
globular clusters remain quite large (±0.3 — 0.9 dex). This reduces the significance of metallicity 
as a linking parameter due to the increased limits of each uncertainty ellipsoid in this dimension. 
Furthermore, it has been demonstrated that the dwarf galaxies of the Local Group exhibit a wide 
variety of star formation histories and metallicities (see Mateo 1998, and references therein). For 
example, Grebel & Guhathakurta (1999) report a large (~ 1 dex) spread in stellar metallicities 
within the M31 satellites And VI and And VII. Similar spreads have been found in the globular 
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cluster systems of some dwarf elliptical galaxies (J. Lotz, private communication). As a result, 
proximity in metallicity is not necessarily a very suitable criterion for group identification within a 
hierarchical formation scenario. 

For these reasons, the complete available sample of 321 M31 globular cluster positions and 
velocities was also examined for grouping without the inclusion of [Fe/H] information. Again, 
roughly a dozen potential groups were identified and then tested for stability against small changes 
in the grouping criteria. The majority of the resulting groups correspond to those found when 
metallicity was included as a linking parameter. 

Perrett et al. (2002) identified significant bimodality in their sample of M31 GC metallicities, 
in confirmation of earlier reports made by Huchra, Brodie & Kent (1991) and Ashman & Bird 
(1993). The centrally concentrated, rapidly rotating, metal-rich cluster system is consistent with 
a bulge population (Perrett et al. 2002). The group-finding results in the full GCS data sample 
both with and without metallicities revealed a couple of relatively large groups that correspond 
to a sizable fraction of this metal-rich bulge population. A third grouping analysis was therefore 
also performed on the sample of 231 GCs in the metal-poor population as defined by Perrett et 
al. (2002). Metallicity was not used as a linking parameter in this test. Several tentative groups 
persisted within the central region even within this metal-poor population, yet the memberships of 
all but one of these groups were not robust to small variations in the linking criteria. Ultimately, 
groups with average galactocentric radii within the central 10' of the galaxy were rejected due to the 
higher probability of false detections, an effect which is confirmed by the Monte-Carlo simulations 
presented later in Section 3.3. Additionally, any cannibalized dwarf GCSs would be extremely 
difficult to recover within the dense inner regions of the galaxy; this effect will also be noted in the 
numerical simulation results discussed in Section 4. 



3.2. Grouping Results 

In the M31 GCS, there is evidence for the presence of ten unique globular cluster groups of 
n GC > 4 members with mean galactocentric radii beyond 10'. The positions of these groups are 
shown in Figure 1, and the data for the group members are provided in Table 1. Five of these 
groups were identified in all three search cases described above. Three additional groups (Groups 
7, 8 and 10) resulted from a search of the full sample of 321 globular clusters without the inclusion 
of [Fe/H] information. Two further groups can be considered as tentative identifications: Group 3 
was found only when metallicity was included as a linking parameter, and Group 4 resulted only 
from the search of the metal-poor population. All of the groups listed in Table 1 are stable to small 
perturbations in the linking criteria. 

An eleventh group is also listed in Table 1; despite the fact that this group lies close to the 
galaxy center, the mean radial velocity of its group members (—155 kms -1 ) is considerably higher 
than the systemic velocity of M31 (—300 kms -1 ). This group deviates significantly from the GCS 
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rotation curve determined by Perrett et al. (2002), and in fact is moving counter to the general 
rotation at that major axis position. 

The association of a group with NGC 205 is clearly evident in Figure 1 (Group 5 of Table 1). 
The mean velocity of the NGC 205 group is v r = —259 kms -1 with a v = 78 kms -1 , a result that 
is consistent with NGC 205's systemic velocity of —242 ± 3 kms" 1 (de Vaucouleurs et al. 1991). 

It is interesting to note that none of the groups recovered in this analysis overlap with the 
associations found by Ashman & Bird (1993). This may be partly attributable to the larger 
sample of higher precision velocities used here. Ashman & Bird (1993) also indicated that their 
method may produce false hits in the presence of GCS rotation. In an effort to reduce this affect, 
they applied a KMM mixture-modelling technique to remove the rotating metal-rich (disk/bulge) 
population from their grouping analysis, thus decreasing their sample from 144 GCs to 95 "halo" 
objects (Huchra, Brodie & Kent 1991). However, Perrett et al. (2002) found that even the metal- 
poor cluster population demonstrates significant rotation, a factor which may also have biased the 
Ashman & Bird (1993) grouping results. 

The majority of the clusters in the inner and outer subsystems found by Saito & lye (2000) 
are included in the metal-rich (bulge) population inferred from the KMM results of (Perrett et al. 
2002). The bulk of inner subsystem found by Saito & lye (2000) was recovered here as two central 
groups; the memberships of these inner groups were not stable and these are not listed in Table 1 
for reasons discussed previously. 

3.3. Monte-Carlo Simulations 

A series of Monte-Carlo simulations were performed to investigate the likelihood that the iden- 
tified M31 globular cluster groups are mere chance associations on the sky. Randomized populations 
of 321 globular clusters were generated based on the globular cluster positions and radial velocities 
from Perrett et al. (2002). Metallicities were not included in the randomized tests due to the large 
[Fe/H] uncertainties in the data. Furthermore, the eight non-tentative outer groups in our sample 
were found without incorporating any metallicity information in the group finding. 

First, the cumulative distribution observed in projected position along the major axis of the 
galaxy was used to select a random cluster X position with roughly the same underlying distribution 
as that seen in the M31 GCS. Next, the cumulative distribution in projected minor axis location of 
the 70 nearest objects in X was used to generate a randomized (but realistic) value of Y position 
for each target. Similarly, the cumulative distribution in radial velocity for the 70 nearest objects 
in X was used to generate a randomized v r for each target. This procedure generated a population 
of globular clusters which shares similar observed properties with the observed data — including 
spatial distribution and rotation profile — while eliminating any small-scale correlations that may 
exist in the dataset. 
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Table 1. M31 globular cluster groups 



ID 


GC 


X (') 


Y (') 


v r (kms 1 ) 


[Fe/H] 


1 


B484-S310 


46.69 


-8.31 


-104 ± 12 


-1.95 ± 0.59 




B480 


44.36 


-8.18 


-135 ± 12 


-1.86 ± 0.66 




DA084 


42.76 


-11.08 


-192 ± 12 


-1.79 ± 0.72 




B376-S309 


42.22 


-10.64 


-142 ± 20 


-2.18 ± 0.99 




B374-S306 


41.14 


-10.50 


-96 ± 12 


-1.90 ± 0.67 


2 


B229-S282 


30.25 


-2.30 


-103 ± 88 


-1.81 ± 0.74 




DAO66 


28.68 


2.88 


-148 ± 12 


-1.82 ± 0.26 




DA065 


27.37 


2.32 


-130 ± 12 


-1.80 ± 0.36 




B216-S267 


26.90 


1.02 


-84 ± 12 


-1.87 ± 0.39 




B223-S278 


26.45 


-3.68 


-101 ± 12 


-1.13 ± 0.51 


3+ 


B203-S252 


21.19 


-0.36 


-199 ± 12 


-0.90 ± 0.32 




B190-S241 


20.95 


2.38 


-86 ± 12 


-1.03 ± 0.09 




B198-S249 


19.98 


0.02 


-105 ± 12 


-1.13 ± 0.30 




B206-S257 


19.75 


-2.23 


-151 ± 23 


-1.45 ± 0.10 




B200 


18.07 


-1.57 


-153 ± 12 


-0.91 ± 0.61 


4+ 


B141-S197 


14.83 


8.06 


-180 ± 12 


-1.59 ± 0.21 




B137-S195 


13.79 


8.50 


-215 ± 12 


-1.21 ± 0.29 




DA058 


13.13 


6.19 


-125 ± 12 


-0.87 ± 0.07 




B102 


12.63 


13.34 


-236 ± 12 


-1.57 ± 0.10 




B105-S166 


9.71 


10.83 


-238 ± 12 


-1.13 ± 0.32 


5 


B11-S63 


3.12 


33.74 


-157 ± 52 


-1.54 ± 0.34 




B328-S54 


3.11 


35.50 


-243 ± 23 


-1.51 ± 0.28 




B324-S51 


2.95 


36.42 


-299 ± 35 


-1.39 ± 0.40 




B9-S61 


1.11 


32.47 


-335 ± 52 


-1.57 ± 0.26 


6* 


B126-S184 


-2.76 


-2.04 


-182 ± 14 


-1.20 ± 0.47 




NB91 


-2.93 


-1.19 


-187 ± 10 


-0.71 ± 0.33 




NB83 


-4.23 


0.84 


-150 ± 14 


-1.26 ± 0.16 




B106-S168 


-4.54 


-0.43 


-90 ± 43 


-0.86 ± 0.68 




B86-S148 


-4.64 


2.48 


-168 ± 21 


-1.74 ± 0.17 


7 


B13-S65 


-7.21 


24.32 


-409 ± 12 


-1.01 ± 0.49 




B16-S66 


-8.99 


21.35 


-406 ± 12 


-0.78 ± 0.19 




B12-S64 


-10.79 


22.98 


-358 ± 12 


-1.65 ± 0.19 




B4-S50 


-11.72 


25.70 


-373 ± 12 


-1.26 ± 0.59 


8 


V31 


-19.06 


7.11 


-433 ± 12 


-1.59 ± 0.06 




V216 


-20.17 


0.96 


-465 ± 12 


-1.15 ± 0.26 




B33-S95 


-21.57 


1.78 


-439 ± 12 


-1.33 ± 0.24 




B31-S92 


-23.12 


1.88 


-400 ± 12 


-1.22 ± 0.40 




B28-S88 


-23.64 


2.54 


-434 ± 12 


-1.87 ± 0.29 




B453 


-23.69 


5.64 


-446 ± 12 


-2.09 ± 0.53 




B15 


-26.57 


7.78 


-460 ± 12 


-0.35 ± 0.96 




DA039 


-26.73 


5.92 


-478 ± 12 


-1.22 ± 0.41 


9 


B57-S118 


-24.94 


-7.16 


-437 ± 12 


-2.12 ± 0.32 




B34-S96 


-26.43 


-2.40 


-540 ± 12 


-1.01 ± 0.22 




B458 


-26.44 


-6.37 


-521 ± 12 


-1.18 ± 0.67 




B49-S112 


-27.49 


-7.41 


-481 ± 12 


-2.14 ± 0.55 




DA048 


-27.91 


-6.55 


-490 ± 12 


-2.01 ± 0.99 


10 


B81-S142 


-25.22 


-12.36 


-430 ± 12 


-1.74 ± 0.40 




B66-S128 


-29.50 


-13.18 


-389 ± 12 


-2.10 ± 0.35 




B65-S126 


-33.24 


-15.81 


-378 ± 12 


-1.56 ± 0.03 




B43-S106 


-33.58 


-11.38 


-414 ± 12 


-2.42 ± 0.51 




B40-S102 


-35.38 


-11.94 


-463 ± 12 


-0.98 ± 0.48 




B342-S94 


-40.35 


-12.27 


-479 ± 12 


-1.62 ± 0.02 


11 


B18-S71 


-40.63 


-4.15 


-585 ± 12 


-1.63 ± 0.77 




B448 


-43.13 


-3.01 


-552 ± 12 


-2.16 ± 0.19 




B335 


-43.95 


-4.87 


-514 ± 12 


-1.05 ± 0.26 




BoD195 


-47.14 


-4.38 


-552 ± 12 


-1.64 ± 0.19 




B327-S53 


-47.65 


-3.51 


-528 ± 12 


-2.33 ± 0.49 




B443 


-50.41 


-4.85 


-532 ± 12 


-2.37 ± 0.46 




B319-S44 


-51.99 


-1.81 


-535 ± 12 


-2.27 ± 0.47 




B315-S38 


-55.62 


-1.13 


-559 ± 12 


-1.88 ± 0.52 





* Central group 
t Tentative group 
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One thousand realizations of the randomized GC populations were performed and each was 
tested for grouping in position and velocity parameter space. The results of these Monte Carlo tests 
are shown in Figure 2. The top panel of Figure 2 shows the radial distribution of the groups found 
in the randomized data. Roughly half of these "false hits" were recovered within the central 10' 
of the galaxy, a region which appears to be too dense to yield reliable group detections. The fact 
that the spatially compact (metal-rich) GC population demonstrates rapid rotation also adds to 
the likelihood of detecting chance conglomerations in position- velocity parameter space. Therefore, 
neglecting all of the groups with mean radii within 10' of the galaxy center, we obtain the number 
distribution of random-data groups shown in the lower panel of Figure 2. 

If we exclude the two tentative groups and the central group, the remaining eight groups 
observed in the actual M31 dataset (Table 1) lie above ~ 97% of the randomized runs. This 
suggests that the substructure found in the M31 GCS is significant overall. If we further exclude 
the group associated with NGC 205, this result still lies above ~ 91% of the randomized runs. 
Therefore, it is unlikely that all of the groups are chance associations, although the majority may 
indeed be. For example, there is nearly a 60% chance that three or more of the groups are not 
real. With the currently available data, it is not possible to determine if the individual groups are 
dynamically bound structures, nor can we gauge the specific degree of internal contamination by 
background clusters. 

4. Modelling the Fate of the Globular Cluster Systems 

Based on the discovery of apparent substructure in the M31 GCS, it is worthwhile to investigate 
the plausibility that the proposed groups are a result of merger events. To accomplish this, it is 
instructive to consider the fate of a dwarf satellite's small globular cluster population during a 
merger with a massive galaxy. N-body simulations were performed to follow the evolution of a dwarf 
galaxy GCS over timescales of r > 1 Gyr as it interacts with an external potential representing 
a massive spiral galaxy. The form of the galactic potential and the modelling of the dwarfs are 
discussed in Sections 4.1 and 4.2, respectively. The interpretation of the simulation results is 
presented in Section 4.3. 

4.1. The Galactic Potential 

In this study, the galactic potential was assumed to be rigid and time- independent, and was 
built from three components: a spheroidal bulge potential (Hernquist 1990), a Miyamoto & Nagai 
(1975) thick disk, and a logarithmic halo. The forms of these potentials were based on the modelling 
of the Sagittarius dwarf by Johnston, Spergel Sz Hernquist (1995) and Helmi & White (2001). The 
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potentials are: 



-GM disk 

<Msk = ; (3) 



^R 2 + (a + Vz 2 + b 2 ) 2 



0halo = Whalo ln ( r2 + ^ 2 ) ( 5 ) 

where G is the gravitational constant, r and represent radii in spherical and cylindrical coordi- 
nates, respectively, and z is the height above the galactic plane. The masses (M), scale lengths 
(a, b, c, d) and velocity (fhalo) m the above potentials were selected to reproduce the characteris- 
tics of the Milky Way. The components of our Galaxy are better understood than those of M31 
and we can expect the overall effects of satellite disruption within M31 to be quite similar. The 
adopted component parameters are provided in Table 2 and have been selected following the work 
of Johnston (1998). 

The representation of the massive galaxy by a fixed potential neglects any exchange of energy 
between the satellite and its parent galaxy. However, typical masses of dwarf satellites are several 
orders of magnitude smaller than that of either the Milky Way or M31. As a result, the effects 
of dynamical friction and energy exchange between the satellite and the massive galaxy can be 
neglected over the timescales of interest here (Helmi & White 2001; Binney & Tremaine 1987). 



4.2. The Model Dwarf 

The dwarf satellite galaxy in each simulation was modelled as a Hernquist sphere with a density 
distribution given by 

M sat s 1 

P{r) = ^rVTs? (6) 



Table 2: Galaxy model parameters 



Parameter 


Adopted Value 


^bulgc 


3.4 x 1O 1O M 


M disk 


1 x 10 n M Q 


^halo 


128 kms" 1 


a 


6.5 kpc 


b 


0.26 kpc 


c 


0.7 kpc 


d 


12.0 kpc 



- 11 - 



(Hernquist 1990), which corresponds to a cumulative mass distribution of 

M{r) = Msat^^. (7) 

M sat is the total satellite mass and s is a scale length that defines the half-mass radius of the 
distribution: setting M(ri/ 2 ) = M sat /2 in Equation 7 gives r 1 / 2 = (1 + Vtys. This analytical mass 
model is an approximation of the de Vaucouleurs i? 1 / 4 profile applicable to galaxies and bulges (de 
Vaucouleurs 1953). 

For each satellite, roughly 10 000 particles were distributed according to this density profile. 
Dwarf galaxies of several masses were considered within the range from 1O 7 M to 1O 9 M . Ten 
of the particles were given a somewhat more extended distribution about the dwarf and were 
designated to be its globular clusters; these GCs were each assigned masses of 5 x 1O 5 M . The 
remaining mass of the simulated dwarfs was equally divided among the rest of the particles. The 
number of particles used to represent the dwarf was sufficient to model its disruption, while being 
small enough to efficiently run the required number of simulations. 

A series of 85 simulations were run in order to explore a reasonable range of parameter space by 
altering the initial conditions of the incident satellite. Table 3 lists the ranges of initial conditions 
that were tested. The total dwarf masses (M sa t) were selected to span a range between a typical 
dwarf spheroidal and a dwarf elliptical and include the mass of its 10 globular clusters. The range 
in satellite half-mass radii (-R1/2) was chosen to correspond to typical values of dwarf galaxy core 
radii — a reasonable assumption if light follows mass. Rgc is the scale radius of the dwarf's 
globular cluster system. The satellites were assigned a starting galactocentric distance (D) and a 
position angle (6) above the plane of the disk. The azimuthal angle was taken to be = 0, since 
the potential is axisymmetric. The parameters a and (3 dictate the initial partition of velocity in 
the XZ and XY plane, respectively, in units of circular velocity such that a 2 + (3 2 = 1 for circular 
orbits. 

These N-body simulations were performed using the Barnes-Hut treecode 1 and modified to 
include an external potential. The treecode algorithm is based on the commonly-used hierarchical 
force method described by Barnes & Hut (1986). 

It should be noted that the simulations presented herein are intended to be illustrative, and 
not to reproduce the detailed characteristics of this complex interaction. The addition of evolving 
potentials for the bulge, disk and halo components and the incorporation of effects such as dynamical 
friction, disk-shocking, and gas dynamics would permit a more comprehensive investigation of the 
merger event itself. 



1 Treecode version 1.4 distributed by J. Barnes. 
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4.3. Grouping in the Simulation Results 

Using the scaling and linking lengths determined in the previous step (A = 0.1, 7„ = 0.13 
km _1 s), the N-body simulation results were analysed to investigate the likelihood of identifying 
the cannibalized dwarf's globular clusters as a group in parameter space. The three-dimensional 
positions and velocities were projected into an orientation similar to M31 (i = 77.7°), and were 
superimposed on the background components described in Section 2.2. The group-finding algorithm 
was then applied to the simulation results at various timesteps. 

Two parameters, Q and Qmax-, were defined in order to quantify the retrieval efficiency of the 
satellite's GCS: 

Qmax = max ( ) (9) 
V rig J k 

where N is the total number of groups containing at least one of the dwarf's original globular 
clusters, uqc is the number of the globular clusters in group k, and n g is the total number of 
objects (including background) that were identified in group k. The sum in Q is calculated only 
over those groups that contain one or more globular clusters and reflects the general contamination 
of the GC groups by background objects. Qmax is the maximum ratio of the number of satellite 
GCs in a group to the total number of objects allocated to that group, and therefore indicates 
whether a bona fide group has been found. For perfect group recovery with no contamination, 

Q = Qmax = 1- 

To demonstrate, let us examine a particular simulation with an initial dwarf position angle 
90° (a = 0.9, = 0) with M sat = 1O 9 M at a timestep of 1.2 Gyr. The group-finding algorithm 
returned globular cluster group memberships as follows: 6 of the 10 dSph GCs were in the same 
group and had 1 additional background object included, 1 GC was allocated to a background group 
with 3 additional members, 1 GC was allocated to a background group with 4 additional members, 



Table 3: Initial conditions for N-body simulations 



Parameter 


Initial Conditions 


M sat 


10 7 - 1O 9 M 


Rl/2 


0.2 - 0.4 kpc 


Rgc 


(1 - 2)i? 1/2 


D 


20 - 100 kpc 


9 


0-90° 


a 


- 0.9 


a 


0-0.9 
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and 2 GCs were not identified as members of any group. For this case with N = 3 groups containing 
one or more globular clusters: 

/6 1 1\ 
Q max = max I -, -, - 1 = 0.86 

The high Qmax indicates that it is very likely that we can detect a group that contains a significant 
fraction of the disrupted satellite's GCS, here with only 14% contamination by background objects. 
The moderate to low Q ratio reflects the fact that some of the satellite's GCs were allocated to 
other groups to which they do not truly belong. 

The orbital periods of the satellite debris within the potential of the parent galaxy are of order 
1 Gyr. The probability of detection clearly depends to a large degree on the relic's location in its 
orbit, as there is a decreased likelihood of finding a group if the GCS debris is fore or aft of the 
bulge, or obscured in the rotating disk particles. Therefore, the likelihood that the detection of a 
group in a galaxy's GCS reveals the remnant of a real accretion event is comparatively low if the 
group is superimposed over the dense part of the disk or the bulge of the galaxy. Not surprisingly, 
satellite GCS remnants are also most easily observable at large galactocentric radii due to their 
significantly lower degree of disruption over the relevant timescales. 

In reality, we rarely know the orbital parameters of past (or even present) satellite galaxies. 
We can therefore examine the average Q values of a large sample of the simulations having a variety 
of initial position angles and incident velocities to get some idea of the likelihood of group detection 
for a satellite on a random orbit. The mean and median Q and Qmax values for satellites with an 
initial galactocentric radius of 40 kpc are shown in Figure 3. The average and median Q values 
decrease with time, as expected. The Q results indicate that it is not likely that we can recover 
the complete GCS of the satellite dwarf in a single uncontaminated group. However, the Qmax plot 
shows that one significant grouping of clusters associated with the cannibalized dwarf can generally 
be found up to several Gyr after the initial encounter. 



5. Conclusions 

The numerical simulation results presented in this study reveal that it is possible to identify 
some fraction of a disrupted satellite's initial GCS based on grouping in projected position and 
radial velocity. This can be accomplished for timescales of r > 1 Gyr after the initial encounter, 
although generally not without some significant degree of background contamination in the group 
membership. 

In the Milky Way, there are three-dimensional data available that trace a distinct stellar 
stream associated with the Sagittarius dwarf merger. It is then relatively easy to pick out any 
globular clusters that coincide with this stream in both distance and velocity, although even here 
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the results are somewhat controversial (for example, see Minniti, Meylan & Kissler-Patig 1996). 
At the distance of M31, we have projected data and are only now beginning to uncover signs of 
stellar streams and substructure (Ibata et al. 2001). These observed streams extend well beyond 
the spatial limits of the current cluster sample (on a line connecting NGC 205 and M32), and thus 
cannot be directly matched with any of the GC groups identified in this study. Large-area surveys 
(e.g., Massey 2002; Lee et al. 2002) will help to provide the necessary additional stellar astrometry 
and to reveal any new GC candidates in M31. Nonetheless, Lynden-Bell &i Lynden-Bell (1995) 
demonstrated that it is possible to use streams of globular clusters to trace the orbits of potential 
merger remnants within the Milky Way. Streams or groups of globular clusters in nearby massive 
galaxies may also provide a viable means of identifying the orbits of any recently merged satellites. 

We have found significant signs of grouping in projected position and radial velocity within 
the globular cluster system of M31. Do these M31 globular cluster groups represent the signatures 
of past accretion events? Without accurate distances or dynamics, the correlation of M31 GCS 
groups with past accretion events remains tenuous. In followup studies, the identification of a 
stellar stream coincident with a grouping of globular clusters would lend the necessary credibility 
to the notion that a galaxy's GCS — and by extension the galaxy itself — formed at least partly 
in a hierarchical manner. 
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Fig. 1. — The positions of the groups identified in the M31 globular cluster system. Globular 
cluster group members are shown with matching symbols and their data are presented in Table 1. 
M31's dwarf companions (NGC 205 and M32) are shown by the large ellipses. 



-18- 



0.3 




0.3 



m 

K 0.2 



c 
o 

3 o.i 

El, 







20 40 60 80 

Average Radius of Group Members (arcmin) 



~i 1 1 1 1 1 1 r 



J_ 



~\ 1 1 1 1 1 r 

> 

GO 



2 4 6 8 10 

Number of Groups Found (r > 10 arcmin) 



Fig. 2. — The grouping results in 1000 randomized populations. The top panel shows the distribu- 
tion in mean galactocentric radius of the false groups found in the randomized data. Roughly 50% 
of the false hits originated within the central 10' of the galaxy (shaded bins). The lower panel shows 
the number distribution of groups with galactocentric radii greater than 10' that were identified in 
the randomized data. The number of groups found in the actual M31 GCS data is shown by the 
dashed vertical line. 
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Fig. 3. — The average Q ratio (top plot) and Q max values (bottom plot) of a series of simulations 
with a range of initial conditions at an initial distance of 40 kpc. The dashed and solid lines show 
smoothed fits to the mean and median values in timestep bins, respectively. 



